Ground-state & finite-temperature properties of spin liquid phase in the J\ - J 2 honeycomb model 
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In this paper we analyze the groundstate and finite-temperature properties of a frustrated Heisenberg J t - J 2 
model on a honeycomb lattice by employing the Schwinger boson technique. The phase diagram and spin gap 
as functions of J2/J1 are presented, showing that the exotic spin liquid phase lies in 0.21 < J2/J1 < 0.43. The 
temperature and magnetic-field dependences of specific heat, magnetic susceptibility and Knight shift are also 
presented. We find the spin liquid state is robust with respect to external magnetic field. These results provide 
clear information characterizing unusual properties of the exotic spin liquid phase for further experiments. 



I. INTRODUCTION 

Frustrated quantum antiferromagnetic systems in low- 
dimensional lattices have received great attention in recent 
years, especially on a honeycomb lattice. Since they have a 
small coordination number (z = 3) in two dimension, large 
quantum fluctuations are expected to lead to novel magnetic 
behaviors. Some magnetic materials with honeycomb lattice 
were synthesized in recent years and exhibit unusual proper- 
ties. A recent high-field electron spin resonance (ESR) spec- 
troscopy measurement on novel Bi3Mn40i2(N03), which is 
a model substance of S = 3/2 honeycomb lattice antiferro- 
magnet (AFM), suggested an important role of the geometric 
frustration in removing the long-range magnetic order, since 
no long-range order was observed down to 1 .9K l . In addi- 
tion, a recent synthesized honeycomb compound In3Cu2VC>9 
also exists weak intralayer and strong interlayer magnetic 
frustrations 2,3 , in which the spin- 1/2 Cu sites forming a honey- 
comb lattice do not exhibit a long-range magnetic order down 
to 2 K. Furthermore, in the family of compounds BaM2(X04)2 
with M=Co, Ni and X=P, Regnault and Rossat-Mignod found 
that the magnetic ions M with spin 5=1/2 for Co or S = 1 
for Ni form frustrated honeycomb layers. These layers sepa- 
rate so large from each other that BaM2(X04)2 can be viewed 
as a two-dimensional AFM 4,5 . The properties of these low- 
dimensional AFMs arise a strong conjecture that whether they 
form a long-searched exotic spin liquid state. 

Frustrated quantum magnets on honeycomb lattice may 
also exhibit remarkable properties. As shown by Meng et al. 6 , 
a quantum spin liquid, which is a short-range resonating 
valence-bond (RVB) liquid due to the frustration induced by 
the motion of the carriers, may be stable in an intermediate- 
interacting Hubbard model on the honeycomb lattice, though 
their results were questioned by Sore 11a et al? . Meanwhile, 
most of theoretical analysis focus on strongly interacting limit 
of the Hubbard model, i.e. the Mott insulating region de- 
scribed by the J\ - J2 AFM Heisenberg model on the honey- 
comb lattice. Soon afterward Meng et al. 's work, Clark et al. 
showed that in a quantum J\ - Jo model a spin liquid phase 
can be stable in the range of 0.08 < J 2 /J\ < 0.3 s . Although 
the variational Monte Carlo approach based on the variational 
family of entangled-plaquette states also supported that a spin 
liquid ground state is stable in the range 0.2 < /2/-/1 < 



0.4 9 , many other authors suggested a plaquette valence bond 
phase in the intermediate /2/-/1 ratio by the renormalization 
group method 10 , the coupled-cluster method 11 and the ex- 
act diagonalizationi 2 - - — . The divergency of zero-temperature 
magnetic phase diagrams on honeycomb lattice obtained by 
different methods shows the necessarity of further study to 
clarify these controversial results. 



Meanwhile, very few work has come to the finite- 
temperature properties of the frustrated honeycomb lattice 
though some have been studied in other structures 15-22 . Such 
theoretical results could provide the experimentists with some 
theoretical guidance to search for exotic spin liquid phase 
in realistic materials. For this purpose, we employ the 
Schwinger-boson mean-field theory (SBMFT), which has 
proved successful in incorporating quantum fluctuation, to in- 
vestigate the finite-temperature properties on the frustrated 
honeycomb lattice. With respect to linear spin-wave the- 
ory, SBMFT can describe both ordered and disordered phases 
even with large fluctuation. Motivated by previous results, in 
this paper we focus solely on antiferromagnetic interactions 
Jij > and present the thermodynamic properties by using 
SBMFT of the Ji - Ji Heisenberg model for spin 5=1/2 



We firstly find that with the increase of J2/J1, the system 
undergoes a transition from a Neel AFM phase to a quantum 
disordered phase at /2/-/1 = 0.21; the disordered phase is clas- 
sified as a spin liquid state, and sequently the disorder phase 
enters into a collinear striped AFM (SAFM) phase at the crit- 
ical point J2/J1 = 0.43 at T - OK. We also clearly show 
that the spin liquid state is robust with respect to strong ex- 
ternal magnetic field; and its low-temperature specific heat, 
magnetic susceptibility and Knight shift exhibit exponential 
behavior. The rest of this paper is organized as follows: we 
first describe the model Hamiltonian and Schwinger boson 
method in Sec. II; for the theoretical and numerical results, we 
present the groundstate properties and zero-temperature phase 
diagram of the Ji - Jo Heisenberg model in Sec. Ill and the 
finite-temperature properties of spin liquid phase in Sec. IV, 
respectively, and the final section is devoted to conclusion. 
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II. MODEL HAMILTONIAN & METHOD 

A. Neel AFM and Disordered Phase 

The honeycomb lattice is composed of two interlacing tri- 
angular lattices. Every site has three nearest neighbors on the 
other lattice, and six next-nearest neighbors on the same lat- 
tice. The J i - J 2 model Hamiltonian reads 

h = J\ y \ s r ■ s r+a + J2 / \ Sr- Sr+p- (l) 

r,a Rfi 

Here J\ and J 2 are the nearest-neighbor and next-nearest- 
neighbor AFM couplings, respectively. This system is frus- 
trated when both J\ and J 2 are positive. The a and [3 sums 
run over all three nearest-neighbor vectors a of site r and the 
six next-nearest-neighbor vectors B of site R. The vectors de- 
scribing the neighbor positions are shown in Fig. Q] Let the 
lattice constant equal to unity. 




FIG. 1: Neel AFM honeycomb lattice with sublattices A and B. The 
vectors a and /? point to the nearest neighbors and next-nearest neigh- 
bors of a central site, respectively. 

Since the honeycomb lattice is a complex Bravais lat- 
tice, we separately treat the sublattice spins with different 
Schwinger bosons. The spin operators are represented as fol- 
lows: 

s t = a l^ a iU S I = «a fl <'T> ( 2 ) 



for sublattice A. The spins on the sublattice B are rotated by n 
from those on the sublattice B. 

a n -* -bju a a -> b jv (3) 

Within the Schwinger boson representation, the sublattice in- 
teractions are rewritten as 

H^—Jx^-.A+jAij-. + lNJiS 2 , (4) 

<ij> 

H 2 = ^J 2 J];B+B il :-^NJ 2 S 2 , (5) 

H 3 = \h J] ■ Cjfijk ■ ~ \nJiS\ (6) 

<j,k> 

where Hi is the Hamiltonian of the interaction between the 
two different sublattices, H2 is that of the sublattice A, and 
H3 is that of the sublattice B; the operators Ay = a^bn + 
a iL bji, Bu = a+ati + ajjjay, and C# = b^b^ + b^bki corre- 
spond to AFM correlation of nearest neighbors and ferromag- 
netic (FM) correlation of next-nearest neighbors in the sublat- 
tice A and B, respectively. N is the number of all lattice sites. 
The total Hamiltonian then reads: 

H -Hi + H 2 + H3 + ^ Ai (a^a^ + aja,-j, - 2S) 
i 

+ Yt*-j( b n b n + b l b *- 2S \ ^ 
j 

where the local constraints for the spin Hilbert space are en- 
forced with the Lagrange multipliers Aj and Aj, which will be 
replaced by a local parameter A in our mean-field treatment. 
Assuming the averages of order parameters of the boson fields 
P =< Ajj > and Q =< Bu >=< Cjk >, we decouple the 
Hamiltonian into quadratic, and the effective H can be diago- 
nalized. 

We use a complex Bogoliubov transformation in the mo- 
mentum space 

aka- = u* k akcr + Vkfitr, bkcr = VjfcOr^ + ulBkcr (8) 

to diagonalize the Halmiltonian by choosing Uk and Vk prop- 
erly. The diagonalized Hamiltonian 

H = J] [Ojkiat^ka + PtJka + l)] + E c (9) 
for 

has the energy dispersion relation of the Schwinger boson ex- 
citation a>k 

= A' - Y l ' (io) 

and the constant energy 

T = \ hpl ~ \ hQl + \ JxSl ~ 3-7252 ~ A ~ 2SA 

and 

X k = A + 3J 2 Qr 2k , Y k = ~7i \Pr lk \ (12) 
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where the geometrical structure factors are given by 



yu 



Thus minimizing the free energy F = -fc^rin Tr(e with 
respect to P, Q and A, we arrive a set of self-consistent equa- 
tions 



-Y 



-(2n* + l) 



(25 + 1) = 0, 



N 



Z 



Iruin 



0>* 



(2m + 1) 



(2n* + 1) 



-P = 0, 

-e = o, 



(14) 
(15) 
(16) 



with n* = l/[exp(jSwj;) - 1] being the Bose-Einstein distribu- 
tion function and /3 = 1/ksT. Solving the mean-field equa- 
tions yields the average values of P, Q, and A, we can explore 
properties of this system at zero and finite temperatures. 

In a two-dimensional honeycomb system, an AFM long- 
range order (LRO) can be stable when Jj is small enough at 
low temperature. As the AFM LRO corresponds to a con- 
densation of the Schwinger bosons, we could obtain the self- 
consistent equations at T — K by the replacement 



S*6(k), 



(17) 



where S * measures the fraction of the Bose-Einstein conden- 
sation (BEC) of the Schwinger bosons 18 . When this system 
is a LRO corresponding to BEC at T = K, the ground-state 
energy per site reads: 



2 v-i E c 

— — > <±>k H 

N ^ N 



(18) 



On the other hand, the system does not undergo the BEC due 
to frustration effect when J 2 /J\ becomes large. We could sep- 
arate the BEC phase with AFM LRO from the non-BEC phase 
with spin disorder by judging the vanishing S * and gap for 
S = 1/2 18 . The ground-state energy per site of the non-BEC 
phase, corresponding to the spin disorder phase, is the same 
to Eq. QD. 



B. Collinear Striped Antiferromagnetic Phase 

When J2/J1 is large enough, from the classical energy anal- 
ysis, one finds that the spin disorder phase becomes unsta- 
ble, and the collinear SAFM phase becomes more stable. The 
magnetic structure and the vectors describing the neighbor po- 
sitions are shown in Fig. [2] 

The SAFM honeycomb lattice composes of four interact- 
ing sublattices (A\, A2, B\ and B2), corresponding to four 
Schwinger bosons. After a rotation of n in spin space in A2 
and Z?2 sites with respect to that in A\ and B\ sites, similar 




FIG. 2: SAFM sublattices A\,A2,B\ and Bi in honeycomb lattice. 
The vectors a , B and 7 point to the nearest and next-nearest neigh- 
bors. 



to Eq. (0, the spin operators in A\ and A2 sublattices can be 
represented as follows: 

s n = 5 ( a h a w ~ a tu a iu) . 

s ti = a m a iu> s i\ = a tu ail i> ( 19 ^ 

S i2 = -a a ^cii2i, S a = -at^a-ai- (20) 

Replacing operator a with b gives rise to the Schwinger 
representation for B\ and B2 sublattices. Rewriting H and 
making the mean-field approximation similar to that in the 
Neel AFM phase, we get an effective Hamiltonian: 



ff = XKl<r Km a k,2<r hoc) 



ka- 

+ E C 

with the constant energy 



(A C* D* B 

C A B D 

D B A C 

B D* C* A) 



"k,2<T 

h + 



(21) 



Ec = 1 
N 2 



+ -J\S + J2S -A- 2SA 
2 



(22) 
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where we define 



A = J 2 Q 2 r lk + A, B = -2J 2 P 2 r 2k , 
C = -JiQir 3k , D = -J\P\r Ak \ 



(23) 



and the geometrical structure factors are 



ri k =e ika \ r 4k 



IZ 



ika\ 



(24) 



Pu Pi, Q\ and Q 2 are the order parameters, corresponding 
to the nearest-neighbor AFM coupling, next-nearest-neighbor 
AFM coupling, nearest-neighbor FM coupling, and next- 
nearest-neighbor FM coupling, respectively, and A is the local 
Lagrange multiplier.After a complex 4x4 Bogoliubov trans- 
formation to diagonalize the Halmiltonian 23 , we get the fol- 
lowing self-consistent equations by minimizing the free en- 
ergy F = -k B T In Tr (e^ H ): 



N 



Z 



»M + 



1 \ dtoi- 



2 dA 



— 

dA 



(25 + 1) = 0, 



(25) 



-y 



n kl + 



N 



Z 



N 



1 \ d(x} k \ I 1 \ du> k2 
2/ <9P! + p 2+ 2/ 3Pi 



+ /iPi = 0, (26) 



l\da) kl I l\do) k2 



+ 2J 2 P 2 = 0, 

(27) 



4 V 17 l\du ki I l\da) k2 ] 1 



(28) 



^ + 2j^ + r + 2)5e7 



-/ 2 62 =0, (29) 



with the Bose-Einstein distribution function n k i, 2 
l/texpfjSwii^) - 1] and the dispersion relations 



= -^A 2 - fi 2 + |C| 2 - \D\ 2 ± y]4A 2 \C\ 2 + 4B 2 \D\ 2 - 4 (ABC*D + ABCD*) + (CD* - CD) 2 



(30) 



Since this LRO state is characterized by a BEC of bosons at 
ko = (0, 0) in the Schwinger boson approach, a straightfor- 
ward calculation yields the ground-state energy per site in the 
collinear SAFM phase 



= ^Z ( ^ i+ ^ 2)+ §-- 



(31) 



Thus, we could compare the ground-state energies among 
Neel AFM, spin liquid (disorder) and SAFM phase, and pick 
out the stablest ground state. 



C. Spin-Spin Correlation Functions 

We could also understand the spin structure of the dif- 
ferent phases by exploring the spin-spin correlation func- 
tions (SSCF). In the presence of the long-range Neel AFM 
order, the SSCF consists of the transverse and longitudinal 
correlations 17 . The longitudinal correlation function between 



different sublattice sites is expressed as 



-2(1 - Mi) 2 - -M\ -(1 -Mi)M 2 , 



where 
M, 



— V — 

N^2co k 



Mi = § Z {[ cos ■ *y) ( r i* + n *) 

+i • sin (^k ■ R*A (r\ k - r lk ) 



3JiP 

8<D k 



(32) 



(33) 



(34) 



The average to various paired terms with four operators in Eq . 

(l32l is in accordance with Wick theorem. 

The transverse correlation functions between sublattice sites 

read 

(5, v • 5*) = (5j • 5-) = -2M 2 + 2M X M 2 - M\ (35) 
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Then in the presence of LRO the SSCFs between the sites 
belonging to the same lattice are calculated, 



(SI ■ Sf) = 2(1 - M) 2 + (1 - Nt)N 2 + \n\ 
(Sf -Sj) = 2N 2 -2N X +N 2 2 , 



where 



2 v— 1 X k 



2tJ k ' 



(36) 
(37) 

(38) 
(39) 



In a similar way we can get the SSCF in the spin liquid phase, 
where the three components of the correlation functions are 
all equal. The SSCF between different sublattice sites thus 
reads 

(^•^•) = -5{|-Z[ cos ( t -^)(^ + '-«) 

+i ■ sin ■ R*^j (r\ k - r lk ) 



3JiP 

8(L> k 



(40) 



The SSCF between different sites belonging to the same sub- 
lattice can be calculated as 



3/ifhti l 

2 



4u> k 

3hP\r lk . 



1 

Aoj u + 2 



(41) 



The SSCF of the collinear SAFM phase can also be obtained 
through calculating longitudinal and transverse correlation 
functions by using the method similar to that of Neel AFM. 



D. Magnetic Field Effect & Knight Shift 

When a magnetic field B is applied parallel to the z-axis, 
similar to the previous process, we can get the diagonalized 
Hamiltonian, 

H M f = ^ (a^a^T + B t[Pk{) + <^kl (a^u +^ T At) 
k 

+2- 



where 



X t = A + 3J 2 Qr 2k , Y k = -h \Pr lk \ 



+ E C (42) 

(43) 
(44) 
(45) 



li = g[*B, g is the Lande factor, hb is the Bohr magneton and 
the constant enegy 



= \ JlPl ~ \ hQl + \ JlSl ~ 3JlSl ~ A ~ 2SA ' (46) 

We get the self-consistent equations by minimizing the free 
energy: 



-7 



("ii + n k2 + 1) 

WkO 



-(25 + 1) = 0, 



N V 

k 

N ^ 



\n k \Y, 



rzkX k 



(«« +n k2 + 1) 



WjfcO 



(n kl +n k2 + l) 



P = 0, 



Q = 0, 



(47) 
(48) 
(49) 



which are reasonable generalization of Eq. (11411161 . Herenti^ 
represents the bosonic occupation with the spectrum 0)41,2- 

Nuclear magnetic resonance(NMR) is a classical method 
for probing the magnetic correlation effects, and Knight shift 
can directly measure the spin fluctuations around a nucleus 
in NMR experiment 24 - 25 . When the magnetic field B parallel 
to z-axis is not too large, the system remains in the spin liq- 
uid phase, the nuclear spin is polarized along this direction, 
and the energy level is split. Simultaneously, the electron spin 
is also polarized, and its polarization can impact on the nu- 
cleus through the hyperfine interaction. Accordingly, NMR 
frequency is changed. The effective magnetic field generated 
by polarized electrons around nucleus is given as follow a 26 ' 27 , 



SB = 



(50) 



where F (q) is the hyperfine interaction structure factor, and 
jm is the nuclear gyromagnetic ratio. Thus, we can get the 
Knight-shift 2 ^, 



K 



Soj 7n-SB F(0) M(T) 



(51) 



y N ■ B y N hfi B 
from which the local magnetic field could be characterized 



III. GROUND-STATE PROPERTIES AND PHASE 
DIAGRAM 

We first present the zero-temperature phase diagram of 
Jz/Ji by comparing the groundstate energy, as shown in Fig. 
[3] There are two critical points. The first one at (J 2 /Ji) c i = 
0.21, which is in agreement with Mattsson's result 2 -, cor- 
responds to a continuous phase transition between the Neel 
AFM and the spin disordered phases. The second one at 
(JilJi)ci = 0.43 corresponds to a first-order phase transition 
between the spin disordered phase and the collinear SAFM 
one. In comparison with the results by Mezzacapo et al., the 
range of our spin disordered phase (0.21 < J 2 /J\ < 0.43) is al- 
most in agreement with that obtained from variational Monte 
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FIG. 4: (Color online) The dispersion relation for J = 0.10 (a), 0.25 
(b), 0.60 (c) and (d) at zero temperature, (c) and (d) correspond to 
two spin waves in the collinear SAFM phase. In the following fig- 
ures, J = JijJ\. 



FIG. 3: (Color online) Ground-state energy per site and energy gap 
in the bosonic dispersion as a function of J2/J1 f° r S = 1/2. For 
J%IJ\ < 0.21 and Ji/Ji > 0.43, the system remains gapless, cor- 
responding to Neel AFM and collinear SAFM phase, respectively. 
When 0.21 < J%IJ\ < 0.43 the system remains gapped. 



Carlo calculations (about 0.2 < J2/J1 < 0.4), based on the 
variational family of entangled-plaquette states 9 . 

In accordance with the phase diagram, the energy gap of 
the Schwinger bosons as a function of J2/J1 is finite only in 
the spin disordered region 0.21 < J2/J1 < 0.43, and the gap 
lifts with the increase of J2/J1, shown in Fig. [3] The energy 
gaps vanish when J2/J1 < 0.21 or > 0.43, where the BEC and 
magnetic LRO appear. In Fig. |4] the spin dispersion relations 
are displayed for Jz/Jy = 0.10, 0.25 and 0.60, respectively. 
It clearly indicates that the BEC occurs at Icq = (0, 0) when 
J2/J1 = 0.10 and 0.60. For J2/J1 = 0.25, the system lies in 
the spin liquid phase. Although the lowest point of the dis- 
persion relation is also at Icq = (0, 0), the first excited state 
is separated from the ground state with a finite value. Hence 
there is no BEC and the system is disordered with a finite gap. 

The SSCF as the function of space distance of spins could 
disclose the spatial distribution of spin arrangement. For 
hfj\ < 0.21, the SSCF is Neel AFM coupling in each direc- 
tion, and for J2/J1 > 0.43 the SAFM long-range correlations 
are also obtained in the zigzag direction, while the SSCF only 



9 
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FIG. 5: (Color online) Spin-spin correlation function as a function of 
distance R in zigzag direction. For the Neel AFM phase J^Ux =0 
and 0.1; for 0.21 < J 2 /J\ < 0.43, the SSCF is short range, indicating 
a gap zone with a short-range order; and for J 2 /J\ > 0.43, the cor- 
relations correspond to the collinear SAFM phase (J = 0.6(A) and 
J = 0.6(6) for two different zigzag directions). 
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shows short-range Neel AFM order for 0.21 < J2/J1 < 0.43 
in the quantum disordered state with a gap in the bosonic 
dispersion. A plot of the SSCF for J 2 /Ji = 0.0, 0.1, 0.25 and 
0.6 in the zigzag direction is shown in Fig. Although Tao 
et al. pointed out that the quantitative features of the present 
SBMFT give the SSCF 1.5 times larger than that of the exact 
result 30 , it is qualitatively correct. For J2/J\ = 0.0, our results 
divided by 1 .5 are in exact agreement with the exact diagonal- 
ization and coupled-cluster method results given by Farnell 14 . 
For /2/-/1 = 0.0 and 0.1, the SSCFs are isotropic and shows 
a classical Neel AFM order behavior. For J%jJ\ = 0.6, the 
SSCFs are different along two different zigzag directions 
due to the SAFM structure, and the long-range correla- 
tions are obtained. When JijJ\ = 0.25, the SSCF rapidly 
decays to zero. Simultaneously, we notice that there is 
no symmetry breaking, so this state is a spin liquid state, 
which can be viewed as a superposition of short-range va- 
lence bond, also called a resonating valence bond (RVB) state. 




FIG. 6: (Color online) Temperature dependence of the RVB order 
parameter P, Q and Lagrane multiplier A for J1/J2 = 4.0. Inset: 
energy gaps with different J1/J2 = 4.0 and 2.5. 



IV. FINITE-TEMPERATURE PROPERTIES OF SPIN 
LIQUID PHASE 

At finite T, the Neel AFM and the SAFM phase become 
disordered phase due to strong gapless spin fluctuations, as 
constrained by the Mermin- Wagner theorem for the present 
two-dimensional system. While the gapped spin liquid phase 
could survive at finite T until up to a critical temperature. 
Thus in the following we are especially interested in the finite- 
temperature properties of the spin liquid phase on the present 
spin frustrated honeycomb lattice. The nearest-neighbor spins 
are paired to form spin singlets, which are destroyed by the 
rising temperature gradually. Therefore, we treat J\, corre- 
sponding to the strength of singlet, as a variable, and let J2 
equal to unity at finite temperatures. 



This results in the linear T behavior of the energy gap when 
the system enters in the paramagnetic phase. 



A. Order Parameters & Energy Gap 

The variations of the order parameter P, Q and Lagrange 
multiplier A with T are obtained by solving the mean-field 
equations for J t /J2 - 4.0. The numerical results are shown in 
Fig. [6] With the increase of the temperature, the spin singlet 
correlations are suppressed by thermal fluctuations gradually. 
When the temperature is increased to about UbTs/Jj = 3.7, 
both of the order parameters P and Q become zero simultane- 
ously. This shows that strong thermal fluctuations in T > T s 
break the spin singlet pairs, and the system transits from the 
quantum disorder phase to the paramagnetic phase. 

Temperature dependence of the energy gap is also presented 
in the inset of Fig. [6] The energy gap, which mainly arises 
from the contribution of the Lagrange multiplier, almost de- 
velops linearly with increasing temperature. Since in the finite 
and high temperatures, the Lagrange multiplier is an effective 
"chemical potential" of the Schwinger bosons; it is consider- 
ably larger than the energy scales of the ordered energies J\P 
and J2Q, and almost linearly increases with the temperature. 



1.6 




'0.0 0.5 1.0 1.5 2.0 2.5 3.0 3.5 4.0 

ksT/Ji 

FIG. 7: ((Color online) Temperature dependence of specific heat per 
site with different /1//2 ' n the spin liquid phase. 
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B. Specific Heat 



In finite temperatures, we calculate the specific heat of the 
spin liquid phase. The shape of C - T is nearly like a A-type 
transition, which is presented in Fig. [7] 

When k B T IJ2 - 0, the specific heat is zero since there is no 
any quasiparticle excitation. In the low-temperature region, 
an analytic form of the low-temperature specific heat per site 
is obtained 



c v 



S C A 2 



Nk B T 



(52) 



where S c is the area of the unit cell of a sublattice, the energy 
gap is 



A = X (3J 2 Q+A) 2 



and the constant 77 is 



9JfP z -36J 2 Q(3J 2 Q+A) 
\6^{3J 2 Q+A) 2 -{\j l pf 



(53) 



(54) 



Here P, Q and A are also constants when T approaches to zero 
temperature. Such a low temperature behavior of spin liquid 
on a honeycomb lattice resembles that on a triangular lattice 16 . 

When T lifts up to T s , the increasing temperature leads to 
more quasiparticle excitations, the contribution of the quasi- 
particle excitations to the specific heat becomes larger and 
larger. Meanwhile, the rising temperature destroys spin sin- 
glet. When T is high enough to T s , the spin singlets van- 
ish and the phase transition from spin liquid to paramagnetic 
phase occurs, as shown in Fig. [7] Different J2/J1 corresponds 
to different strengths of spin singlet, leading to different criti- 
cal temperatures. When the system transits from a RVB spin 
liquid state to a paramagnetic state, the specific heat sharply 
drops to zero, corresponding to a second-order phase transi- 
tion similar to that seen in a square lattice^. 



C. Magnetic Susceptibility 

Under different magnetic fields, we get the temperature de- 
pendence of the RVB order parameter and Lagrange multi- 
plier, which is presented in Fig. |8]for J\ /J 2 - 4.0. Through- 
out this paper, we keep external magnetic field IJ.BIJ2 < 1.39, 
actually, which is not a small value. When = 1.39, we 

obtain B = 60 T by putting the value of J2 = 5 meV. We find 
that the influence of the magnetic field on the order parame- 
ters is rather small, suggesting the spin liquid phase is robust 
and self -protected to applied magnetic field, as seen in Fig. [8] 
This arises from that the spin singlet in the RVB spin liquid 
phase is rigid with respect to external magnetic field B until a 
critical pair-breaking field B c . 

In finite magnetic fields, the energy dispersion of the 
Schwinger bosons splits into two branches, hence two energy 
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FIG. 8: Temperature dependence of the order parameters P, Q and 
Lagrange multiplier A in the magnetic field B = and 3, respectively. 
In the following figures, B is measured in unit of fiB/Jo. 



gaps. The magnetic field dependence of the energy gaps, cor- 
responding to Eq. (|43l at k = (0,0), is shown in Fig. [9] for 
J1/J2 = 4.0 and k B T/J 2 = 0.345. With the increase of the 
magnetic field, two energy gaps go in different ways, one gap 
lifts and the other decreases. The separation between the two 
gaps linearly increases with applied magnetic fields, about en- 
ergy scale of fiB. 

We also investigate the temperature dependence of magne- 
tization and susceptibility under different magnetic fields for 
J1/J2 - 4.0. The low-temperature behavior of the magnetic 
susceptibility per site is given as follows: 



x(T )^ill e -Nk B T CQsh MB 



87T77 



2k B T 7 



(55) 



where the energy gap A and constant 77 are expressed in Eq. 
( 153b and ( f54b . respectively. When B approaches to 0, we get 
the magnetic susceptibility analytically 



X ~ 8^ 6 



(56) 



These low-T behaviors on a honeycomb lattice are very simi- 
lar to Mila et oZ.'s on a square lattice 18 . There are three tem- 
perature ranges in Fig. [10] Firstly, in extremely low tem- 
perature an exponential behavior of magnetic susceptibility as 



9 



presence of the short range coupling and effective magnetic 
field lead to a Curie-Weiss law. 



0.0 0.2 0.4 0.6 0.8 1.0 12 1.4 

/JS/J2 



FIG. 9: (Color online) The separation of two energy gaps at k B T / Ji 
0.345 for J x lh = 4.0. 



expressed in Eq. ( pol l. In the second temperature range, the 
magnetic susceptibility almost increases linearly with the lift 
of the temperature. In < T < T s , the spin singlet formed 
between two neighboring spins is gradually destroyed by the 
thermal fluctuations. Hence the spin is easier to be polarized 
so that the magnetization increases with the lift of tempera- 
ture. Since the magnetic susceptibility characterizes the de- 
gree of difficulty for polarization and magnetization, it also 
becomes large for higher temperature when T < T s . Finally, 
for T > T s , the system transits from the spin liquid phase to 
the paramagnetic one. With the increase of the temperature, 
the thermal fluctuations of spins become more and more in- 
tense, depressing the spins polarization. Therefore, the mag- 
netization and magnetic susceptibility gradually decrease in 
the high temperature region. Moreover, the inverse suscepti- 
bility is linearly dependent on temperature and consistent with 
the Curie law. Similar behavior is also obtained by Mila et al. 
for square lattice 31 , demonstrating that the SBMFT also works 
well in high temperature T > T s . 

One may notice that such a paramagnetic phase with the 
Curie law is completely different from the conventional para- 
magnetic phase with the Curie- Weiss law. The quantum spin 
fluctuations in the former leads to vanishing static magnetic 
field on central spin, hence to the Curie law in magnetic sus- 
ceptibility, suggesting a quantum paramagnetic phase. For 
the latter a conventional paramagnet with AFM coupling, the 




FIG. 10: (Color online) Temperature dependence of magnetization 
(a) and susceptibility (b) per site in different magnetic field B for 
J1/J2 = 4.0. Inset in (b): the reciprocal of the susceptibility (c) in 
different magnetic field B for J[IJi = 4.0. 



D. Knight Shift 

The low-temperature behavior of the Knight shift per site is 
given by 



K 

K — — 

N 



F(0) S c k B T c _ 
y^ti 4jttjB 



A/k B T 



sinh 



(57) 



where the energy gap A and constant rj are expressed in Eq. 
d53l ) and d54l >, respectively. When B approaches to 0, its ana- 
lytic form reads 



F(0) _ A^ e _, 



(58) 



Namely the Knight shift is exponentially small at T — > 0. 

The temperature dependence of the Knight shift for differ- 
ent Ji/J 2 at yuB/7 2 = 0.46 is shown in Fig. [TT] With the 
increase of temperature, the Knight shift lifts with tempera- 
ture with different slope for different ratio J\/J2 in the spin 
liquid phases. It is interestingly found that the smaller the ra- 
tio J1/J2 is, the larger the Knight shift is in the spin liquid 
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FIG. 11: (Color online) Temperature dependence of Knight shift 
per site for different J1IJ2 at /jB/Ji = 0.46. The unit of K is 
F(Q)/(y N frfj) . Inset: temperature dependence of Knight shift per 
site in different magnetic field for J1JJ2 = 4.0. 



phase. This is because the strength of spin singlet is weak for 
small J\/J2, and the electron spin can be easily polarized by 
applied magnetic field, leading to a large Knight shift. Once 
the spin liquid-paramagnetic phase transition occurs, all of the 
T-dependent Knight shifts fall into a single curve in the param- 
agnetic phase. These behaviors could be easily understood as 
follows: for T < T s , the rising temperature destroys the spin 
singlet and the AFM spin correlations become weak gradually. 
Accordingly, the electron spin is easier to be polarized so that 
the effective magnetic field around nuclear becomes larger. 
Thus the nuclear magnetic resonance (NMR) frequency and 
the Knight shift lift with the increasing temperature in T < T s . 
When T > T s , the system enters the paramagnetic phase, the 
intense thermal fluctuation considerably reduces the effective 
magnetic field around a nucleus. Hence, the Knight shift de- 
creases with the increase of temperature in the paramagnetic 
phase. 

We also present the magnetic field influence on the Knight 
shift of spin liquid phase. The temperature dependence of 
Knight-shift for different magnetic field at J1/J2 = 4.0 is 
shown in the inset of Fig. QT| We find that the Knight-shift 
almost does not change when fiB varies from zero to 1 .39/2- 
This is attributed to the fact that the effect of magnetic field on 
the order parameters P and Q is small, as seen in Fig. [8] which 
leads to the spin singlet almost unchanged. Thus the Knight- 
shift in different magnetic fields changes very little. From Fig. 
QT|and the inset one notices that the critical temperatures ob- 
tained from the Knight shift are also in agreement with those 



obtained from the specific heat in Fig. [7] 



V. REMARKS 

From the preceding study we have found that through com- 
paring the total energy, there are three stable phases in the 
phase diagram of frustrated J\ - J2 Heisenberg model on a 
honeycomb lattice: the Neel AFM phase in J2/J1 < 0.21, 
the RVB spin liquid phase in 0.21 < J2/J1 < 0.43, and 
the collinear SAFM phase in J%jJ\ > 0.43 . Although sev- 
eral authors suggested a plaquette valence bond crystal state 
in the intermediate frustration regime, instead of the RVB 
spin liquid one, through the exact diagonalization method^— 
and coupled-cluster method^, it was not known whether the 
finite-size result is stable in the thermodynamic limit ( i.e. 
N ~ infty ). We also notice that the critical value of the spin 
liquid-Neel AFM boundary in the Hubbard model 8 is smaller 
than that of ours, this is attributed to the duplicate effects of 
the spin fluctuations and the charge fluctuations, which pre- 
vent the formation of magnetic long-range order. 

We find that beyond the phase transition temperature T s , 
the magnetic susceptibility recovers the usual Curie behavior. 
These results demonstrated that the Schwinger-boson mean- 
field theory works well for calculating thermodynamic prop- 
erties of the frustrated J\ — J2 Heisenberg model on a hon- 
eycomb lattice, even if in the high temperature regime with- 
out magnetic order. This arises from the finite value of the 
gap parameters A, even if the order parameters P and Q van- 
ish in the high temperature. More importantly, we present 
the low-temperature exponential behavior of the specific heat, 
magnetic susceptibility and Knight shift, providing a clue for 
searching the candidate compounds with exotic spin liquid 
phase in experiments. We also illustrate that it is not easy for 
an external magnetic field destroy the spin liquid phase, show- 
ing a self-protecting and robust characteristics of this quantum 
phase. 
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